{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d7fae1bb-fbe6-4eb0-8faf-de4d26a781ad",
   "metadata": {},
   "source": [
    "# Survival Analysis TSR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81745670-55df-4a9e-b961-87c6eb248c9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import csv\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import train_test_split, KFold, cross_val_score\n",
    "from sklearn import metrics\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n",
    "from sklearn.metrics import auc\n",
    "from sklearn.metrics import RocCurveDisplay\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "s\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47405e79-fff7-47e5-9628-e2bc95399f04",
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(x,axis=None):\n",
    "    x_max = np.amax(x, axis=axis, keepdims=True)\n",
    "    exp_x_shifted = np.exp(x - x_max)\n",
    "    return exp_x_shifted / np.sum(exp_x_shifted, axis=axis, keepdims=True)\n",
    "    \n",
    "def display_auc(X1,X10,X100,y1,y10,y100,name=''):\n",
    "    cv = KFold(n_splits=5, shuffle=False)\n",
    "\n",
    "    tprs_old = []\n",
    "    aucs_old = []\n",
    "    old_predictions, new_predictions,avg_predictions = [],[],[]\n",
    "    old_labels,new_labels,avg_labels = [],[],[]\n",
    "    tprs_new = []\n",
    "    aucs_new = []\n",
    "    \n",
    "    tprs_avg = []\n",
    "    aucs_avg = []\n",
    "    mean_fpr = np.linspace(0, 1, 100)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(15,10))\n",
    "\n",
    "    for i, (train, test) in enumerate(cv.split(X1, y1)):\n",
    "        classifier_old = LogisticRegression()\n",
    "        classifier_old.fit(X1[train], y1[train])\n",
    "        y_pred = classifier_old.decision_function(X1[test])\n",
    "        my_prediction = classifier_old.predict_proba(X1[test])\n",
    "        my_prediction = softmax(my_prediction, axis=1)[:,1]\n",
    "        old_predictions.append(my_prediction)\n",
    "        old_labels.append(y1[test])\n",
    "        fpr, tpr, thresholds = metrics.roc_curve(y1[test],y_pred)\n",
    "        roc_auc = metrics.auc(fpr, tpr)\n",
    "        viz = metrics.RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=roc_auc,\n",
    "                                      estimator_name=\"ROC fold {}\".format(i))\n",
    "\n",
    "        interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)\n",
    "        interp_tpr[0] = 0.0\n",
    "        tprs_old.append(interp_tpr)\n",
    "        aucs_old.append(viz.roc_auc)\n",
    "\n",
    "\n",
    "    for i, (train, test) in enumerate(cv.split(X10, y10)):\n",
    "        classifier_new = LogisticRegression()\n",
    "        classifier_new.fit(X10[train], y10[train])\n",
    "        y_pred = classifier_new.decision_function(X10[test])\n",
    "        my_prediction = classifier_new.predict_proba(X10[test])\n",
    "        my_prediction = softmax(my_prediction, axis=1)[:,1]\n",
    "        new_predictions.append(my_prediction)\n",
    "        new_labels.append(y10[test])\n",
    "        fpr, tpr, thresholds = metrics.roc_curve(y10[test],y_pred)\n",
    "        roc_auc = metrics.auc(fpr, tpr)\n",
    "        viz = metrics.RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=roc_auc,\n",
    "                                      estimator_name=\"ROC fold {}\".format(i))\n",
    "\n",
    "        interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)\n",
    "        interp_tpr[0] = 0.0\n",
    "        tprs_new.append(interp_tpr)\n",
    "        aucs_new.append(viz.roc_auc)\n",
    "    \n",
    "    for i, (train, test) in enumerate(cv.split(X100, y100)):\n",
    "        classifier_avg = LogisticRegression()\n",
    "        classifier_avg.fit(X100[train], y100[train])\n",
    "        y_pred = classifier_avg.decision_function(X100[test])\n",
    "        my_prediction = classifier_avg.predict_proba(X100[test])\n",
    "        my_prediction = softmax(my_prediction, axis=1)[:,1]\n",
    "        avg_labels.append(y100[test])\n",
    "        avg_predictions.append(my_prediction)\n",
    "        fpr, tpr, thresholds = metrics.roc_curve(y100[test],y_pred)\n",
    "        roc_auc = metrics.auc(fpr, tpr)\n",
    "        viz = metrics.RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=roc_auc,\n",
    "                                      estimator_name=\"ROC fold {}\".format(i))\n",
    "\n",
    "        interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)\n",
    "        interp_tpr[0] = 0.0\n",
    "        tprs_avg.append(interp_tpr)\n",
    "        aucs_avg.append(viz.roc_auc)\n",
    "\n",
    "    ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=2, color=\"r\", label=\"Chance\", alpha=0.8)\n",
    "\n",
    "    mean_tpr_old = np.mean(tprs_old, axis=0)\n",
    "    mean_tpr_old[-1] = 1.0\n",
    "    mean_auc_old = auc(mean_fpr, mean_tpr_old)\n",
    "    std_auc_old = np.std(aucs_old)\n",
    "    ax.plot(\n",
    "        mean_fpr,\n",
    "        mean_tpr_old,\n",
    "        color=\"b\",\n",
    "        label=r\"TSR 6 months = %0.2f $\\pm$ %0.2f\" % (mean_auc_old, std_auc_old),\n",
    "        lw=2,\n",
    "        alpha=0.8,\n",
    "    )\n",
    "\n",
    "    mean_tpr_new = np.mean(tprs_new, axis=0)\n",
    "    mean_tpr_new[-1] = 1.0\n",
    "    mean_auc_new = auc(mean_fpr, mean_tpr_new)\n",
    "    std_auc_new = np.std(aucs_new)\n",
    "    ax.plot(\n",
    "        mean_fpr,\n",
    "        mean_tpr_new,\n",
    "        color=\"r\",\n",
    "        label=r\"TSR 12 months = %0.2f $\\pm$ %0.2f\" % (mean_auc_new, std_auc_new),\n",
    "        lw=2,\n",
    "        alpha=0.8,\n",
    "    )\n",
    "    mean_tpr_avg = np.mean(tprs_avg, axis=0)\n",
    "    mean_tpr_avg[-1] = 1.0\n",
    "    mean_auc_avg = auc(mean_fpr, mean_tpr_avg)\n",
    "    std_auc_avg = np.std(aucs_new)\n",
    "    ax.plot(\n",
    "        mean_fpr,\n",
    "        mean_tpr_avg,\n",
    "        color=\"y\",\n",
    "        label=r\"TSR 18 months = %0.2f $\\pm$ %0.2f\" % (mean_auc_avg, std_auc_avg),\n",
    "        lw=2,\n",
    "        alpha=0.8,\n",
    "    )\n",
    "\n",
    "    std_tpr_old = np.std(tprs_old, axis=0)\n",
    "    tprs_upper_old = np.minimum(mean_tpr_old + std_tpr_old, 1)\n",
    "    tprs_lower_old = np.maximum(mean_tpr_old - std_tpr_old, 0)\n",
    "    ax.fill_between(\n",
    "        mean_fpr,\n",
    "        tprs_lower_old,\n",
    "        tprs_upper_old,\n",
    "        color=\"blue\",\n",
    "        alpha=0.2,\n",
    "        label=r\"$\\pm$ 1 std. dev.\",\n",
    "    )\n",
    "\n",
    "    std_tpr_new = np.std(tprs_new, axis=0)\n",
    "    tprs_upper_new = np.minimum(mean_tpr_new + std_tpr_new, 1)\n",
    "    tprs_lower_new = np.maximum(mean_tpr_new - std_tpr_new, 0)\n",
    "    ax.fill_between(\n",
    "        mean_fpr,\n",
    "        tprs_lower_new,\n",
    "        tprs_upper_new,\n",
    "        color=\"red\",\n",
    "        alpha=0.2,\n",
    "        label=r\"$\\pm$ 1 std. dev.\",\n",
    "    )\n",
    "    \n",
    "    std_tpr_avg = np.std(tprs_avg, axis=0)\n",
    "    tprs_upper_avg = np.minimum(mean_tpr_avg + std_tpr_avg, 1)\n",
    "    tprs_lower_avg = np.maximum(mean_tpr_avg - std_tpr_avg, 0)\n",
    "    ax.fill_between(\n",
    "        mean_fpr,\n",
    "        tprs_lower_avg,\n",
    "        tprs_upper_avg,\n",
    "        color=\"yellow\",\n",
    "        alpha=0.2,\n",
    "        label=r\"$\\pm$ 1 std. dev.\",\n",
    "    )\n",
    "\n",
    "    ax.set(\n",
    "        xlim=[-0.05, 1.05],\n",
    "        ylim=[-0.05, 1.05],\n",
    "        title=\"Receiver operating characteristic\",\n",
    "    )\n",
    "    ax.legend(loc=\"lower right\")\n",
    "#     plt.savefig('/data/pathology/projects/Pdac-epithelium-segmentation/Survival/Survival_all/shuffled_combined_'+name+'.png')\n",
    "    plt.show()\n",
    "    return np.concatenate(old_predictions),np.concatenate(old_labels),np.concatenate(new_predictions),np.concatenate(new_labels),np.concatenate(avg_predictions),np.concatenate(avg_labels)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0009bef0-5a48-4891-bcf0-d34df0e0f7ee",
   "metadata": {},
   "source": [
    "## Prepare paths and dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b062797-9e35-4793-b52c-cb60f742966b",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical_path = '/data/pathology/projects/Pdac-Segmentation/clinical.cart.2021-12-13/Clean_clinical.tsv'\n",
    "tsr_path = \"/data/pathology/projects/Pdac-epithelium-segmentation/TSR_Final/Tsr_combined_TCGA.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4688a127-16ee-4455-82e8-ca709e0e2db2",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsv_data = pd.read_csv(clinical_path, sep='\\t')\n",
    "df1 = tsv_data[['Patient_id', 'Survival label', 'days_to_survival','gender', 'age_at_index', 'primary_diagnosis','prior_malignancy','origin']]\n",
    "df1 = df1.drop_duplicates(subset='Patient_id', keep='last')\n",
    "tsr_data_old = pd.read_csv(tsr_path, sep=',')\n",
    "tsr_old_names = tsr_data_old['case_id']\n",
    "tsr_old_names = [i[:12] for i in tsr_old_names]\n",
    "tsr_old_names = np.unique(tsr_old_names)\n",
    "\n",
    "assert len(tsr_old_names) == len(np.unique(tsr_old_names)), print('There are {} duplicates'.format(len(tsr_old_names)-len(np.unique(tsr_old_names))))\n",
    "# df2 = tsr_data_old[['case_id', 'tsr_ratio', 'median_tsr', 'average_tsr']]\n",
    "df2 = tsr_data_old[['case_id','tsr_ratio_CH']]\n",
    "df2=df2.rename(columns={\"tsr_ratio_CH\": \"tsr_ratio\"})\n",
    "\n",
    "df2['case_id'] = df2['case_id'].apply(lambda x: x[:12])\n",
    "\n",
    "\n",
    "df_names = pd.DataFrame(tsr_old_names, columns =['case_id'])\n",
    "# df_names = pd.DataFrame(tsr_old_names, columns =['case_id'])\n",
    "df2 =pd.merge(df_names,df2, on='case_id' )\n",
    "df2 = df2.sort_values(by='tsr_ratio')\n",
    "df2 = df2.drop_duplicates('case_id', keep='first')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afd7490-aab9-4327-a7c9-953ca056730c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1=df1.rename(columns={\"Patient_id\": \"case_id\"})\n",
    "df1=df1.rename(columns={\"Survival label\": \"vital_status\"})\n",
    "df3 = df1.drop_duplicates(subset='case_id', keep='last')\n",
    "my_df1 = pd.merge(df3, df2, on='case_id')\n",
    "my_df1.loc[my_df1['vital_status']=='Alive', 'vital_status'] = 0\n",
    "my_df1.loc[my_df1['vital_status']=='Dead', 'vital_status'] = 1\n",
    "my_df1['age_at_index'] = my_df1['age_at_index'].apply(lambda x: x/100)\n",
    "my_df1.loc[my_df1['gender']=='male', 'gender'] = 1\n",
    "my_df1.loc[my_df1['gender']=='female', 'gender'] = 0\n",
    "\n",
    "my_df1.loc[my_df1['prior_malignancy']=='no', 'prior_malignancy']                = 0\n",
    "my_df1.loc[my_df1['prior_malignancy']=='yes', 'prior_malignancy']               = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b6189c-3000-41d4-9eb7-cefc59ecd0b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_df_categorical_old = my_df1.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24d57ff1-618f-436a-b4f9-fc58e5148e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "primary_diagnosis = pd.get_dummies(my_df_categorical_old['primary_diagnosis'], drop_first=True)\n",
    "tissue_or_organ_of_origin = pd.get_dummies(my_df_categorical_old['origin'], drop_first = True)\n",
    "my_df_categorical_old.drop(['primary_diagnosis', 'origin'], axis=1, inplace=True)\n",
    "my_df_categorical_old = pd.concat([my_df_categorical_old, primary_diagnosis, tissue_or_organ_of_origin], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88a63b2b-ab89-41ee-b4d2-32cb3b3ced59",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "my_df_categorical_old.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "682de63a-c325-48a9-b30c-cc5ea6f3f377",
   "metadata": {},
   "source": [
    "## Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b90bf45a-f482-44b6-b1e5-8f84fb854ff6",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold  = 180\n",
    "threshold_2 = 365\n",
    "threshold_3 = 540\n",
    "\n",
    "my_df_categorical_old['state_6'] = my_df_categorical_old.apply(lambda row: 1 if (row['vital_status'] == 1) and (row['days_to_survival'] <= threshold) else 0, axis=1)\n",
    "my_df_categorical_old['state_12'] = my_df_categorical_old.apply(lambda row: 1 if (row['vital_status'] == 1) and (row['days_to_survival'] <= threshold_2) else 0, axis=1)\n",
    "my_df_categorical_old['state_18'] = my_df_categorical_old.apply(lambda row: 1 if (row['vital_status'] == 1) and (row['days_to_survival'] <= threshold_3) else 0, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "388af18e-02da-4219-9c34-493edbb9be55",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_df_categorical_old.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9607377e-258b-4551-a385-91c55f1c2045",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_categorical_shuffled = my_df_categorical_old.sample(frac=1, random_state=42)\n",
    "df_categorical_shuffled.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f5b994-b300-4515-bba7-1faf223c1187",
   "metadata": {},
   "outputs": [],
   "source": [
    "feats = ['gender', 'age_at_index', 'prior_malignancy', 'tsr_ratio', 'Infiltrating duct carcinoma, NOS', 'Pancreas, NOS', 'Tail of pancreas']\n",
    "\n",
    "Xc1 = df_categorical_shuffled[feats]\n",
    "yc1 = df_categorical_shuffled['state_6']\n",
    "Xc1 = np.asarray(Xc1)\n",
    "yc1 = np.asarray(yc1).astype('int')\n",
    "\n",
    "\n",
    "yc10 = df_categorical_shuffled['state_12']\n",
    "yc10 = np.asarray(yc10).astype('int')\n",
    "\n",
    "\n",
    "yc100 = df_categorical_shuffled['state_18']\n",
    "yc100 = np.asarray(yc100).astype('int')\n",
    "\n",
    "\n",
    "short_pred,short_labels,one_y_pred,one_y_labels,eighteen_m_pred,eightee_m_labels =display_auc(Xc1,Xc1,Xc1,yc1,yc10,yc100,name='Features')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e7dc31f-9afd-47e2-b53f-ff530389de0c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
